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I.  INTRODUCTION 

When  the  contact  surface  of  an  expanding  gas  moves  faster  than  the  speed 
of  sound  in  the  stagnant  medium  ahead  of  it,  a  blast  wave  is  formed.  Solutions 
to  this  type  of  problem  are  of  interest  in  a  variety  of  practical  problems 
such  as  the  rupture  of  pressurized  vessels  or  the  blast  from  high  energy 
explosives.  Ignoring  chemical  reactions  within  the  blast  wave,  the  rate  at 
which  work  is  done  in  pushing  aside  the  gas  ahead  of  the  contact  surface  is 
the  prime  factor  in  determining  which  analytical  model  will  most  accurately 
describe  the  real  blast  wave.  For  instance,  a  self-similar  point  source 
model  such  as  that  used  by  Taylor1  assumes  that  all  the  work  of  expansion  is 
done  initially,  in  an  infinitesimal  amount  of  time,  and  is  valid  only  for 
strong  blast  waves  (i.e.,  high  shock  Mach  numbers).  A  more  realistic  self¬ 
similar  model  such  as  that  used  by  Freeman2  assumes  continuous  energy  input 
to  the  blast  wave  but,  except  for  the  case  of  a  constant  wave  Mach  number,  is 
also  valid  only  in  the  strong  blast  regime.  Clearly  what  is  needed  is  an 
analytical  technique  which  is  valid  for  all  blast  wave  strengths  and  uses  the 
realistic  assumption  that  the  energy  input  exists  for  the  finite  amount  of 
time  during  which  the  piston  (contact  surface)  expands.  Such  a  technique 
is  presented  here;  it  is  nonsimilar  and  assumes  an  asymptotic  energy  input 
to  the  blast  wave.  The  method  is  applied  to  the  specific  case  of  the  rupturing 
of  a  sphere  of  compressed  gas  and  the  results  compared  with  the  numerical  -- 
method  of  characteristics  --  solution  to  the  identical  problem. 


1 1 .  PROCEDURE 

Self-similar  solutions  assume  that  the  shock  normalized  blast  wave 
variables  do  not  have  explicit  time  dependence.  Though  this  assumption  may 
be  useful  when  the  wave  is  strong,  it  must  be  inaccurate  for  the  entire  blast 
wave  history.  For  example,  the  minimum  in  the  blast  wave  pressure  will  be 
located  at  its  innermost  boundary  when  the  wave  is  strong  but  moves  forward 
as  the  wave  overexpands.  Self-similar  models  correctly  predict  the  strong 
wave  behavior  but  without  intrinsic  time  dependence  are  unable  to  predict  a 
change  in  the  minimum  location  as  the  wave  weakens.  The  method  used  here  is 
nonsimilar  and  is  a  modification  of  the  procedure  used  by  Bach  and  Lee3 
which  is  not  used  intact  because  of  its  assumption  of  an  instantaneous  point 
source  of  energy.  The  Bach  and  Lee  analysis  is  modified  here  to  treat  the 
more  realistic  case  of  an  asymptotic  energy  supply. 


1.  G.  I.  Taylor ,  "The  Formation  of  a  Blast  Wave  by  a  Very  Intense  Explosion, 
I,  Theoretical  Discussion ,"  Proceedings  of  the  Royal  Society  (London), 
Series  A,  Vol.  201,  1950,  pp.  159-174. 


2. 


R.  A.  Freeman,  "Variable-Energy  Blast  Waves,"  Brit.  J.  Appl.  Phys. 
(J.  Phys.  D),  Ser.  2,  Vol.  1,  1968,  pp.  1697-1710. 


3. 


G.  G.  Bach,  and  J. 
AIAA  Journal,  Vol. 


H.  S.  Lee,  "An  Analytical  Solution  for  Blast  Waves,  " 
8,  No.  2,  February  1970,  pp.  271-275. 


Convenient  coordinates  for  a  spherically  symmetric  blast  are 


y  =  r/Rs  ,  (1) 

where  R  is  the  shock  front  radius,  and 
s 

"  -  (‘A)2  .  C2) 

where  cq  is  the  stagnant  medium  sound  speed,  with  the  dot  symbolizing  time 

differentiation.  Let  the  density,  pressure,  and  fluid  velocity  within  the 
blast  wave  be  defined  by,  respectively, 

p(r,t)  =  pQ  A(y,n)  ,  (3) 

p(r,t)  =  P0(RS)2  B(y,n)  #  (4) 

u(r,t)  =  Rg  C(y,n)  .  (5) 

One  of  two  a  priori  assumptions  associated  with  this  analytical  technique  is 
to  assume  that  the  density  in  the  blast  wave  conforms  to 

P(r,t)  =  p(Rs,t)  yq  ,  (6) 

or  equivalently 

A(y,n)  =  A(l,n)  yq  ,  (7) 

where  q  is  a  function  of  n  to  be  determined.  Assuming  (7)  for  A  makes  it 
possible  to  solve  for  B  and  C  from  the  equations  for  mass  and  momentum 
conservation  alone,  viz. 

(C-y)A  +  ACy  +  2CA/y  =  2hnAn  ,  (8) 

(C-y)C  +  hC  +  B  /A  =  2hnC  ,  (9) 

v  '  y  y'  n 

where 

h  =  Rs'R’s/(^s)2  (10) 

with  the  subscripts  y  and  n  indicating  partial  differentiation  with  respect 
to  y  or  n  accordingly.  The  Rankine-Hugoniot  boundary  relations  across  the 
shock  front  are,  for  a  polytropic  gas, 

Al  =  A(l,n)  =  (y+1)/ (y-l+2n)  ,  (11) 

Bj  =  B(l,n)  =  (2Y  +  n(l -Y) )/ (Y (Y+1) )  ,  (12) 
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Cx  =  C(l,n)  =  2(l-n)/ (Y+l)  , 


(13) 


where  y  is  the  specific  heat  ratio  in  the  blast  wave  region.  Substitution 
of  (7)  for  A  into  the  continuity  Equation  (8)  yields  upon  integration  the 


following  closed  form  expression  for  C: 

C(y,n)  =  y(T  +  (CrT)y"(q+3)  +  Ylny)  ,  (14) 

where 

T(q,n,h)  =  G(q)  +  hH(n,q)  ,  (15) 

Y(q,n,h)  =  hJ(n,q)  ,  (16) 

G(q)  =  (q)/(q+3)  ,  (17) 

H(n,q)  =  (2nA^) / (Aj (q+3) )  +  J/(q+3)  ,  (18) 

J(n,q)  =  (2nq ')/ (q+3)  ,  (19) 


with  the  prime  signifying  total  differentiation  with  respect  to  n.  Substi¬ 
tution  of  (14)  for  C  and  (7)  for  A  into  the  momentum  Equation  (9)  yields 
upon  integration  a  similarly  closed  form  expression  for  B, 

B(y,n)  =  Bj  +  Ax  (2hnT'-hT-(T-l) (T+Y)) (yq+2-l)/(q+2) 

+  A1(2hn(C>1-T")-h(C1-T)+(T-l)(C1-T)(q+2)-(C1-T)(T+Y))(l-l/y) 

+  \  (2hnY '-hY-Y (2T+Y-1) ) ( (lny-1/ (q+2) ) (yq+2/ (q+2) ) +1/ (q+2) 2) 

+  A1(2hn(T-C1)q'-Y(T-C1) (q+1) ) (1- (lny+1) (1/y)) 

+  Aj (Y2) (2/ (q+2) 3- (ln2y- (2/ (q+2) ) (lny-1/ (q+2) ) (yq+2/ (q+2) ) 

+  A1((C1-T)2(q+2))(l-y"(q+4))/(q+4)  .  (20) 


The  right-hand  side  of  the  closed  form  fluid  variable  expressions,  (7),  (14), 
and  (20),  are  completely  known  at  any  given  time  provided  the  quantities 
q (t) ,  h(t),  and  n(t)  are  known  at  that  time.  The  variables  q,  h,  and  n 
are  determined  by  invoking  global  mass  and  energy  conservation  across  the 
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blast  wave  region  as  follows.  Conservation  of  mass  within  the  blast  wave 
requires  that 


j  rv  _ 

/  s  4Ttpr  dr  =/  s  4TTpQr  dr  , 


where  R  and  R  (t)  are  the  contact  surface  radius  at  time  t=0  and  time  t, 
o  c 

respectively.  Figure  1  shows  the  blast  wave  region,  defined  between  R^  and  R^, 

as  it  evolves  in  time.  Substitution  of  (6)  into  (21)  yields  upon  integration 
an  algebraic  equation  for  q: 

q  =  3Aj (l-wq+^)/ (1-wq^)  -  3  ,  (22) 


where 


w  =  R  /R 
c  s 


w  =  R  /R 
o  o  s 


Using  (22)  for  q  introduces  the  unknown  Rc(t);  note  Rq  will  be  given  as  an 
initial  condition  of  the  problem.  To  obtain  R^(t)  a  second  a  priori 

assumption  is  provided  (the  first  being  (6));  viz.,  it  is  assumed  in  advance 

that  an  expression  for  R  (t)  exists.  The  following  expression  for  R  (t)  was 

c  c 

found  useful  for  modeling  the  contact  surface  of  propellant  plumes  from  guns 

and  will  be  used  here  (though  the  theory  remains  valid  for  any  contact 

surface  expression) : 


R  (t)  =  R  (1-lM+t/x)  +  R 
c  00  0 


where  R^  and  x  are  two  constants  which  for  the  present  can  be  assumed  known. 

Note  that  expression  (25)  for  the  contact  surface  radius  is 
asymptotic  in  time.  With  Rc(t)  known,  the  expression  (22)  for  q  depends 

solely  on  n  and  w.  The  three  quantities  h,  n,  and  w  are  all  related 
through  the  time  derivatives  of  Rg(t).  By  definition  w  satisfies 

w  =  w(Rc-wco/v/n) /Rc  .  (26) 

Similarly,  the  inverse  Mach  number  squared,  n,  obeys 

n  =  -2c  hwv^n/R  .  (27) 

0  c 
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Finally,  an  equation  for  h  is  developed  from  global  energy  conservation  as 
follows.  For  a  perfect  gas,  let  the  internal  energy  per  unit  mass  be  given 
by 

e  =  P/(P(Y-1))  .  (28) 

Using  (28),  global  energy  conservation  is  written 

Rs(t) 

E  (t)  =  /  (p/(Y-l)+pu2/2)47Tr2dr 

Rctt) 

-  4TrRj(l-wo3)(poco2/(3Y(Y-l)))  ,  (29) 

where  Es(t)  is  the  work  input  to  the  blast  wave  from  the  pushing  action  of 
the  piston-like  contact  surface  and  is  therefore  given  by 

Rs(t) 

E  (t)  =  /  p(w)47rr2dr  .  (30) 

R 

o 

Integration  and  rearrangement  of  (29)  using  (30)  yields  for  dh/dt 

h  =  A(L-W)/V  ,  (31) 

where 

W(h,n,q,w)  =  E  (t)/((4TTP  c  2R  3)(l-w  3))  (32) 

o  U  U  b  O  9 

L(h,n,q,w)  =  -(l-wo3)/(3k(k-l))  +  ((l-w2)/2) (-a4- (ag/2) +b2- (b5/2) ) 

+  ((1-w3)/3) (a1-a2+a4+a6+ag+2a9+a10)  +  ((l-w'Cq+1))/(q+l)) (a1Q+b4) 

-  ((w(q+5)ln2w)/(q+5))(b6-(q+2)2a9)  +  ((l-w(q+5))/(q+5)) (a2-a6(2q+7)/(q+5) 

-  2a^( (2q+7)/ (q+5) +(q+2) 2/ (q+5) 2)  +  b1-(b3/(q+5))+(2b6/(q+5)2)) 

+  (w2lnw/2)  (ag-b5)  -  (w(-q+5-)lnw/(q+5))  ((q+2)a& 

+2ag  ((q+2)(2q+7)/(q+5))+b3-(2b6/(q+5)))  , 


(33) 


V  (h,n,q,w)  =  a5((l-w  )/2)  +  ((lV)/3)  (a3-a5-a?) 

+  ((l-w(q+5))/(q+5)Ha7(2q+7)/(q+5)-a3)  +  ( (w(q+5)  lnw) / (q+5) )  (a? (q+2) ) ,  (34) 

with  the  a^'s,  b^'s  and  c^'s  defined  in  the  appendix.  The  h  dependence  in 

(32)  for  w  can  be  factored  out  and  moved  to  the  left-hand  side  in  (31)  but 
is  left  undone  here  to  shorten  the  expressions.  It  is  convenient  for 
numerical  solution  to  put  (22)  for  q  into  the  form 

q=n(((-2A1(l-wl-q+3-))/(Y-l  +  2n))-A1(q+3)  (w(-q  +  2^)w/n) )/ (1-wq3) 

+  3wo3A1(l-w(q+3))/((2hn)(l-wo3)2))/((A1w(q+3)lnw)/(l-wo3)  ♦  1/3)  .  (35) 

The  system  of  four  first  order  ordinary  differential  equations  for  q, 
h,  n,  and  w  (viz.  (35),  (31),  (27),  and  (26))  can  be  solved  numerically 
using  a  Runge-Kutta  procedure  provided  initial  values  for  q,  h,  n,  and  w 
are  available.  To  obtain  these  it  is  assumed  that  the  initial  expansion 
of  the  contact  surface  is  analogous  to  the  initial  expansion  after  rupture 
of  a  shock  tube  diaphram.  The  theory  for  the  shock  tube  case  has  a  well- 
known  analytic  solution.  It  will  thus  be  assumed  for  simplicity  that  the 
gas  expands  for  an  initial  time  interval  AtQ  as  if  it  were  a  one-dimensional 

planar  problem.  Values  for  R  (At  ),  R  (At  ),  R  (At  ),  and  R  (At  )  are  thus 
v  r  sv  o’  c  o  s  o’  sv  o 

obtained  and  used  in  (22),  (23),  (2),  and  (10)  to  provide  initial  values  for 

q,  w,  n,  and  h.  Having  obtained  the  solution  for  q(t),  h(t),  and  n(t),  their 

substitution  into  (7),  (14),  and  (20)  yields  the  blast  wave  density,  velocity, 

and  pressure,  respectively,  at  that  time  as  a  function  of  y  alone. 

This  treatment  is  fundamentally  different  from  References  1-3  in  that  it  is 

not  a  point  blast  theory.  Unlike  References  1-3,  R  ft)  and  R  are  not  assumed 

c  o 

zero.  In  fact,  replacing  Rc(t)  and  Rq  by  zero  in  the  above  expressions 

reduces  them  to  their  point  blast  counterparts  in  Reference  3.  The  present 
method  is  now  demonstrated  on  a  hypothetical  test  case. 


III.  RESULTS 


Using  (25)  for  Rc(t),  assume  that  an  ideal  gas,  with 

y  (specific  heat)  =  1.25  , 

T  (temperature)  =  2831  °K  , 

v  (moles)  =  1.63  moles  , 

m  (mass/mole)  =  0.024  kg/mol  , 


E 


« 


T 

!V 

H 

t 


ft 


is  confined  within  a  spherical  volume,  of  radius 

R  =  0.1  m  (37) 

o  » 

which  ruptures  at  time  t  =  0.  The  solution  for  the  ensuing  blast  wave  is 
predicted  using  both  the  present  analytical  technique  and  a  more  exact  non- 
isentropic  method  of  characteristics  technique.  For  comparison,  the  results 
are  also  nondimensionalized  and  compared  with  the  point  blast  model  of  Bach 
and  Lee. 

For  the  purpose  of  evaluating  R  ,  used  in  (25) ,  let  it  be  assumed  that 

the  compressed  gas  undergoes  an  irreversible  and  adiabatic  expansion  against 
a  background  pressure  of  one  atmosphere.  Under  these  conditions,  conservation 
of  energy  will  lead  to  the  following  equation  for  the  final  radius 

R*  =  (3vkT/(4TTYP0))(1/3)  -  Rq  ,  (3 

where  k  is  the  universal  gas  constant  and  pQ  is  one  atmosphere  of  pressur 

both  in  the  SI  system.  Furthermore,  to  evaluate  t,  in  (25),  use  again  the 
assumption  that  the  initial  expansion  of  the  compressed  gas  in  the  radial 
direction  can  be  approximated  by  a  one-dimensional  planar  expansion  in  a 
shock  tube.  In  which  case  the  contact  surface  velocity  at  t  =  0  is  given  by 

Rc(t=o)  =  Rj{2x)  -  CoZ/(Y(1+(Y^1)Z/(2y))1)  ,  (39) 


where  the  first  equality  in  (39)  is  obtained  from  the  time  derivative  of  (25) 
at  t  =  0,  and  Z  in  the  last  equality  must  satisfy 


IS 

Z/(Y(1  +  (Y+1)Z/(2Y))*)  =  (2c/(co(y-l)))(l-(Po(l+Z)/p)a'1)/2Y)  , 

(40) 

• 

where 

« 

p  =  vkr/(4TTR03/3)  , 

(41) 

h 

i 

fc.’ 

c  =  (YP/(vm/(4wRo3/3)))i  . 

(42) 

► 

V 

From  (39) 

* 

T  =  Rco(Y(1+(Y+1)Z/(2y))1/(2coZ)) 

(43) 

1  Though  the  analysis  leading  to  (38)  and  (43)  is  oversimplified,  the 

k  expressions  thus  obtained  were  useful  for  the  intended  purpose  of  (25),  viz., 

•  to  predict  the  contact  surface  radius  of  the  exhaust  gases  expelled  after 

shot  ejection  from  guns.  In  any  case,  the  expression  used  for  the  contact 
’  surface  radius  is  not  at  issue  in  the  present  treatment;  it  is  merely  assumed 


;  ii 

r 

1 1 
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a  priori  that  such  an  expression  exists.  Obtaining  (38)  and  (43)  for  use  in 
(25)  gives  the  appropriate  R£(t)  for  the  example  under  investigation.  The 

complete  results  for  this  test  case  are  shown  in  Figures  2-4.  Figure  2  is  a 

composite  of  the  time  variation  of:  (a)  the  exponent  q  used  in  (6)  for  the 

power  law  form  of  the  density;  (b)  the  Mach  number  of  the  shock  front, 

1 

M  =  n  1 ;  and  (c)  the  fraction  of  energy  input  at  time  t  to  the  total  energy 
input  at  time  t  =  “.(The  scale  of  time  is  chosen  to  show  the  most  dramatic 
change  in  the  variables,  which  in  this  problem  is  from  t  =  0  to  t  =  4500usec.) 
Figure  2(c)  shows  that  the  energy  input  to  the  blast  wave  is  indeed  asymptotic 
with  time  as  opposed  to  the  instantaneous  input  characteristic  of  point  blast 
theories.  A  comparison  between  the  results  predicted  by  the  point  blast 
theory  of  Bach  and  Lee  and  the  present  treatment  is  summarized  in  non- 
dimensional  format  in  Figure  3.  Perhaps  surprising  in  view  of  their  inherent 
differences  is  the  closeness  in  the  curves  of  Figure  3(a)  for  the  shock 
trajectory  R  /R  vs  n>  where 

,  1/3 

Re  *  (Etotal/(«,co  4"»  (44> 

It  appears  from  this  plot  that  approximating  a  timed  release  of  energy  E  ^ 

by  an  instantaneous  release  of  Etotal  has  essentially  no  effect  on  the 

predicted  shock  trajectory.  On  the  other  hand,  Figures  3(b)  -  3(d)  show  that 
there  are  differences  in  the  predicted  profiles  of  the  other  fluid  variables. 
(The  subscript  1  refers  to  the  value  at  the  shock.)  Most  noticeable  of  these 
is  the  difference  in  the  inner  boundary  of  the  blast  wave  which  is  r  =  0  for 
the  point  blast  model  and  r  =  R  (t)  for  the  piston  blast  model.  In  Figures 
3(c)  and  3(d),  the  piston  blast  solution  falls  within  the  limits  predicted 
by  the  point  blast  solution.  But  Figure  3(b)  shows  a  significant  difference 
between  the  piston  and  point  blast  solutions.  In  the  latter  case,  the  fluid 
velocity  has  a  substantial  negative  phase  while  in  the  former  case  the 
velocity  is  only  beginning  to  show  a  region  of  reversed  flow.  In  order  that 
Figure  3  may  be  taken  with  confidence  to  illustrate  the  inherent  limitations 
in  the  point  blast  approximation,  a  comparison  is  made  between  the  piston 
solution  used  in  Figure  3  and  a  more  exact  numerical  solution.  To  serve 
this  purpose,  a  nonisentropic,  numerical  method  of  characteristic  (M.O.C.) 
solution  was  obtained  for  this  problem.  A  comparison  between  solutions  is 
shown  at  four  different  times  in  Figure  4.  The  particular  times  chosen 
illustrate  collectively  the  full  range  of  variable  profiles.  The  time 
spanned  in  Figure  4  correlates  with  the  time  scale  in  Figure  2.  In  addition, 
the  n  =  0.82  solution  of  Figure  3  corresponds  to  the  analytical  solution 
shown  in  Figure  4(c).  Shown  at  each  time  in  Figure  4  are  four  subplots. 

The  time  of  Figure  4(a), for  example, is  indicated  by  the  dashed  horizontal 
line  in  the  top  subplot.  Following  this  dashed  line  across  to  the  piston  and 
shock  curves  and  then  down  to  the  horizontal  axis  (distance  axis)  defines  the 
blast  wave  region  at  this  time.  The  blast  wave  region  is  marked  through  each 
of  the  lower  subplots  by  vertical  dashed  lines.  Kach  of  the  lower  subplots  - 
pressure,  density,  and  velocity  -  shows  the  profile  of  the  corresponding  fluid 
variable  as  it  changes  from  the  piston  to  the  shock.  The  analytical  solution 
shows  good  agreement  in  all  four  subplots  when  the  blast  is  strong,  Figures 
4(a)  and  4(b).  The  favorable  comparison  has  diminished, however,  at  the  time 


of  Figure  4(c).  Figure  4(d)  shows  that  the  disagreement  which  began  in  the 
density  variable,  cf.  Figure  4(c),  now  appears  in  all  the  fluid  variables 
except  the  shock  trajectory,  which  as  noted  in  Figure  3  is  a  relatively 
insensitive  quantity. 


IV.  CONCLUSION 

The  analytical  point  blast  solution  of  Bach  and  Lee,  which  is  based  on 
a  power  law  form  for  the  density,  is  extended  here  to  describe  the  more 
realistic  case  of  blast  waves  driven  by  an  asymptotic  piston-like  motion  of 
an  expanding  inner  gas.  For  the  same  total  blast  energy,  it  was  found  that 
the  piston  and  point  blast  solutions  predict  virtually  the  same  shock  tra¬ 
jectory.  However,  the  point  blast  model  predicts  overexpansion  in  the 
velocity  profile  significantly  earlier  than  the  piston  blast  solution.  The 
analytical  piston  blast  solution  compared  well  with  the  more  exact  M.O.C. 
solution  at  Mach  numbers  greater  than  approximately  1.1.  Theoretical 
disagreement  appears  first  in  the  density  profile  but,  as  the  wave  weakens, 
is  observed  in  all  the  fluid  variables,  with  the  exception  of  the  shock 
trajectory.  Thus,  the  analytic  technique  presented  here  for  blast  waves 
driven  by  a  finite  piston  expansion  supersedes  its  point  blast  predecessor, 
Reference  3,  by  giving  a  more  accurate  description  of  real  blasts.  The 
technique  is  also  a  viable  alternative  to  numerical  solutions  provided  the 
Mach  number  of  the  blast  wave  remains  above  approximately  1.1. 


D  Lstance 

Fig.  i  Blast  wave  pressure  contours  at  two  widely  separated  times,  where 
R  (t)  is  the  contact  surface  radius,  R$(t)  is  the  shock  front  radius,  Rq 

is  the  initial  (preburst)  sphere  radius,  and  PQ  is  the  ambient  background 


pressure 


Fig.  2  Time  variation  of:  a)  the  density's  power  low  exponent  q,  b)  the 
shock  Mach  number,  and  c)  the  fractional  energy  input  to  the  blast  wave. 
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APPENDIX  A 


APPENDIX  A 


The  a.'s,  b.'s  and  c.'s  used  in  (32)  -  (34)  are  as  follows: 

11  l 

=  BjVOnCy-l)) 

a2  =  A1((2hn)(G%hH')-hT+(l-T)(T+Y))/(n(q+2)(Y-l)) 
a3  =  A12hH/((q+2)(y-l)) 

a4  =  A1((2hn)(C'  -G^hH^)fh(T-C1)  +  (T-l)(C1-T)(q+2)  +  (T-C1)(T+Y))/(n(Y-l)) 
a5  =  -A12hH/(Y-l) 

a6  =  A1(2h2nJ'-hY+Y(l-2T-Y))/(n(Y-l)(q+2)2) 
a?  =  A12hJ/((Y-l)(q+2)2) 

ag  =  A1(2hn)(T-C1)(q/n)+Y(C1-T)(q+l))/(n(Y-l)) 
ag  =  A1Y2/(n(Y-l)(q+2)3) 
a10  =  ( (q+2) (Cj-T) 2)/(n( Y-l) (q+4) ) 

bl  =  A1T2/2n 

b2  =  A1T(c1-T)/n 
b3  =  A1TY/n 

b4  =  -A1CC1-T)2/2n 

b5  =  AjY^j-D/n 

b£  =  A.Y2/2n 
o  1 
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